数学模型

Wang Haihua

🍈 🍉🍊 🍋 🍌

---

定积分的计算

定积分的计算是Monte Carlo方法引入计算数学的开端，在实际问题中，许多需要计算多重积分的复杂问题，用Monte Carlo方法一般都能够很有效地予以解决，尽管Monte Carlo方法计算结果的精度不很高，但它能很快提供出一个低精度的模拟结果也是很有价值的。而且，在多重积分中，由于Monte Carlo方法的计算误差与积分重数无关，因此它比常用的均匀网格求积公式要优越。

问题1

求二重积分 $I=\iint_{x^{2}+y^{2} \leq 1} \sqrt{1-x^{2}} d x d y$.

求解

根据积分的几何意义, 它是以 $\sqrt{1-x^{2}}$ 为曲面顶, 以 $z=0, x^{2}+y^{2} \leq 1$ 为底的柱体 $D$ 的体积。用下列简单思路求 $I$ 的近似值, $D$ 被包含在长方体 $\Omega$ : $|x| \leq 1,|y| \leq 1,0 \leq z \leq 1$ 的内部, 长方体 $\Omega$ 的体积为 4 。而 $D$ 是 $x^{2}+y^{2} \leq 1$, $0 \leq z \leq \sqrt{1-x^{2}}$ 所围成的区域。若在 $\Omega$ 内产生均匀分布的 $N$ 个点, 有 $n$ 个点落 在 $D$ 的内部。由频率近似于概率, 得到在 $\Omega$ 内任取一点, 落在 $D$ 内的概率 $p=\frac{D \text { 的体积 }}{\Omega \text { 的体积 }} \frac{I}{4} \approx \frac{n}{N}$, 所以 $I \approx \frac{4 n}{N}$, 计算得 $I \approx 2.6666$ 。

代码

from numpy.random import uniform
import numpy as np
N=10000000; x=uniform(-1,1,size=N)
y=uniform(-1,1,N); z=uniform(0,1,N)
n=np.sum((x**2+y**2<=1) & (z>=0) & (z<=np.sqrt(1-x**2)))
I=n/N*4; print("I的近似值为：",I)

---

from numpy.random import uniform
import numpy as np
N=10000000; x=uniform(-1,1,size=N)
y=uniform(-1,1,N); z=uniform(0,1,N)
n=np.sum((x**2+y**2<=1) & (z>=0) & (z<=np.sqrt(1-x**2)))
I=n/N*4; print("I的近似值为：",I)

I的近似值为： 2.6653972

x1 = np.linspace(-1,1,1000)
x2 = np.linspace(1,-1,1000)
y1 = np.sqrt(1-x1**2)
y2 = -np.sqrt(1-x2**2)

x = np.concatenate([x1.reshape((1,-1)),x2.reshape((1,-1))],axis=1).flatten()
y = np.concatenate([y1.reshape((1,-1)),y2.reshape((1,-1))],axis=1).flatten()

X,Y = np.meshgrid(x,y)
Z = (X**2+Y**2)
import matplotlib.pyplot as plt
flg = plt.figure('zfunc')
ax = flg.add_subplot(projection='3d')
ax.plot_surface(X, Y, Z)
ax.set_xlable = ('x')
ax.set_ylable = ('y')
plt.show()

ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['bottom'].set_position(('data',0))
ax.spines['left'].set_position(('axes',0.5))
plt.plot(x,y)
plt.fill_between(x1,y1,y2,alpha=0.3)
plt.axis('equal')
plt.show()